Coupled online sequential extreme learning machine model with ant colony optimization algorithm for wheat yield prediction

Inadequate agricultural planning compounded by inaccurate predictions results in an inflated local market rate and prompts higher importation of wheat. To tackle this problem, this research has designed two-phase universal machine learning (ML) model to predict wheat yield (Wpred), utilizing 27 agricultural counties’ data within the Agro-ecological zone. The universal model, online sequential extreme learning machines coupled with ant colony optimization (ACO-OSELM) is developed, by incorporating the significant annual yield data lagged at (t − 1) as the model’s predictor to generate future yield at 6 test stations. In the first phase, ACO is adopted to search for suitable, statistically relevant data stations for model training, and the corresponding test station by virtue of a feature selection strategy. An annual wheat yield time-series input dataset is constructed utilizing data from each selected training station (1981–2013) and applied against 6 test stations (with each case modelled with 26 station data as the input) to evaluate the hybrid ACO-OSELM model. The partial autocorrelation function is implemented to deduce statistically significant lagged data, and OSELM is applied to generate Wpred. The two-phase hybrid ACO-OSELM model is tested within the 6 agricultural districts (represented as stations) of Punjab province, Pakistan and the results are benchmarked with extreme learning machine (ELM) and random forest (RF) integrated with ACO (i.e., hybrid ACO-ELM and hybrid ACO-RF models, respectively). The performance of the ACO-OSELM model was proven to be good in comparison to ACO-ELM and ACO-RF models. The hybrid ACO-OSELM model revealed its potential to be implemented as a decision-making system for crop yield prediction in areas where a significant association with the historical agricultural crop is well-established.

The current wheat yield prediction and forecasting methods adopted by the Pakistan Government have been reported to be highly inaccurate 12 . In 2005, poor wheat yield predictions of Pakistan where the actual production was relatively small equated to the estimated yield, resulted in an inflated local market rate and prompted higher importation of wheat 13,14 . Similarly, in 2012-2013, Pakistan experienced severe challenges in wheat supply, which happened due to the lower production of wheat yield in Punjab 15 . A plausible reason for this deficit was attributed to poor agricultural planning and inaccurate predictions to satisfy the national grain needs. Sajjad reported about the looming wheat scarcities for an agriculturally rich country, Pakistan 16 . Due to these concerns that have a direct detrimental impact on income and food security for the already staggering economy of developing Pakistan, the government representatives work towards enhancing the forecast to account for the surplus and shortfalls in advance.
The modelling of wheat yield utilizing ancient procedures and irrelevant information from past yields is bound to deter the outcomes drastically 17,18 . To calculate the future productions, establishing novel systems for improving current and future agricultural productions, and supporting future food security issues in both developing and first-world nations are necessary. This validates the essential role of wheat yield modelling via novel artificial intelligence (AI) modelling or data-intelligent approaches that encapsulate relevant historical patterns. Data-intelligent techniques have enormous flexibility in crop management due to their ease of employment, viable accuracy, and feature detection capabilities [19][20][21][22] . These models are also bound to empower officials in attaining efficient ways to predict future crop productions 23 . There are several examples of data intelligent algorithms in agronomy. Dempewolf et al. implemented vegetation index to predict wheat yield in Punjab 24 , whereas Hamid examined the wheat frugality and future forecasts 25 . On the other hand, Muhammad investigated the historic context of the wheat improvements for Balochistan 26 . Furthermore, Iqbal et al. designed an autoregressive integrated moving average model (ARIMA) to predict future wheat and yields up to the year 2022 in Pakistan 27 . Sher and Ahmad integrated the Cobb-Douglas function with ARIMA to predict wheat yield 28 . However, the works have constructed simplistic regression models (e.g., ARIMA) that are often discredited due to their assumptions of linearity in the data 29 .
Based on previous approaches for wheat yield prediction, Rahman et al. developed a data-driven approach to predict rice yield for Bangladesh 30 , whereas monitoring of rice crop was implemented via a neural network model 31 . Similarly, an artificial neural network (ANN) was developed for soybean and corn predictions in Malaysia 32 . In addition, all the forgoing works were on a provincial level, or a country-wide, which lacks the significance to a small locality such as the district level forecasting used in this study for better accuracy and applicability. Yield prediction is a challenging job as various interconnected climatic drivers affect the yield 5 . Thus, agronomic experts can possibly use the preceding yields to predict future production. Despite this, none of the previous work has utilized wheat yield of several locations for training purposes to predict the yield of other stations.
Information from several other locations for training purposes to predict the yield at the main region is essentially beneficial in decision support systems since it can allow the modellers to accept analogous features prevailing to be analysed to evaluate the main region data 33 . This framework can be adopted in agricultural practices by associating station-specific crop production and creating suitable deductions relevant to the existence of favourable (or unfavourable) eco-friendly or soil fertility circumstances to produce maximum yield 34 . Considering the need for accurate future wheat yield prediction, the modelling of crop yield using several stations yield data for model development can offer a reasonable system to determine the most cost-effective and useful agricultural monitoring practices.
The proposed two-phase AI system called (i.e., ACO-OSELM) model was adopted in the current research to predict wheat yield. Two benchmark models including random forest (RF) and extreme learning machine (ELM) and their hybrid versions were designed for verification of the ACO-OSELM model. ELM and RF models were selected as a benchmark due to their remarkable predictive potentials as appear in the literature [35][36][37][38][39] . The selection of the OSELM was owing to the main merit of the ELM model 40 . ELM model is a single layer feed-forward neural network (SLFN) where the input weights are randomly assigned while the output weights are analytically determined 41 . Unlike conventional neural network models, the ELM is able to avoid issues such as tuning of learning rates, learning epochs, stopping criteria and local optima making it computationally efficient 42 . In addition, ELM efficiently handles large-scale data with a better generalization capability and is more suited to largescale wheat yield predictions 43,44 . On the other hand, the RF models are ensemble regression tree models that use the bootstrap aggregation (i.e., bagging) approach to generate forecasts [45][46][47] . The RF model ameliorates the overfitting issue, which is a key drawback of conventional solitary regression tree-based models 48 . Consequently, these models have been applied in this study. The two-phase hybrid ACO-OSELM is validated for wheat yield prediction in agronomic regions: Rahimyar Khan, Dera Ghazi Khan (denoted as D. G. Khan), Kasur, Sialkot, Rawalpindi, and Jhang located in Punjab province, Pakistan where 26 stations from 26 districts were used to develop the model. The selected study stations are spread throughout the Punjab province and are the major wheat producers. These stations are chosen randomly from the agriculturally rich Punjab province. To verify the applicability of the proposed ACO-OSELM model, this study aims to fulfill three objectives: (i) To develop a bio-inspired ACO algorithm to select the best possible neighbouring stations located in Punjab province, Pakistan for training purposes using feature selection strategy; (ii) To incorporate the statistically important one step antecedent data (i.e., t − 1 where t represents the current data) of the selected training stations in the OSELM model to develop a two-phase hybrid ACO-OSELM model to predict the current and future wheat yield; and (iii) To assess the predictive accuracy of the two-phase ACO-OSELM model for wheat yield prediction universally in the whole province of Punjab in Pakistan.

Theoretical overview
The architecture involved in the establishing of a two-phase hybrid ACO-OSELM model for wheat yield prediction is discussed here.
Ant colony optimization (ACO) algorithm. Dorigo and Di Caro 49 presented the ACO feature selection procedure, which has been widely used in different applications [50][51][52][53] . This study utilized the ACO technique to determine the least possible distance between wheat yield (W) of the training stations, and the testing stations, a strategy that can be adopted to choose the respective training stations for yield prediction at the testing station. A parameter named pheromone in the ACO process is allotted to predictor stations, which categorizes these predictors alongside the target/test stations at the start. The trial pheromone value is used to compute the probability of the chosen station to train against the test station while the magnitude pheromone alters by navigating the training stations, and subsequently, the probability is improved for the next coming ants to pick the optimum station. Readers can survey the following literature for more details on the ACO procedure experiment 54,55 . Extreme learning machine (ELM). Huang  The term G(a, b, x) is representing a nonlinear piecewise continuous function satisfying ELM universal approximation capability theorems 57 , (a, b) are the hidden node parameters and R is the set of real numbers whereas R d is the d-dimensional set of real numbers and x is the input data. The activation functions are Sigmoid, Hyperbolic tangent, Gaussian, Hard limit, Cosine and Fourier basis functions.
Initially, ELM randomly modifies the hidden layer to project the inputs into a feature space using some piecewise continuous nonlinear functions 58 . The parameters (a, b) are generated randomly that are not dependent on the training set. In the second phase of ELM learning, then, the weights (ρ) linking the hidden and the output layer are solved by minimizing the prediction error in the squared error sense: i.e.
where M is denoting the hidden layer output matrix and T is the training data matrix which can be simplified as follows 57 . The ∥ · ∥ indicates the Frobenius norm.
The ideal solution to (3) is provided by: In Eq. (6) M + is indicates the Moore-Penrose generalized inverse of M. The SLFNs with randomly chosen input weights successfully learn various training patterns with the least error. In this way, SLFNs can be considered as a linear system. The output weights which attach the hidden layer to the output layer in the linear system can now be analytically solved by generalized inverse operation of the hidden layer output matrices. Thus, the ELM model is faster than the conventional feedforward learning algorithms 59,60 . Online-sequential extreme learning machine (OSELM). In OSELM the data is channelled in a chunk-by-chunk manner for better understanding and accuracy, whereas in ELM a total of N training data points are used for training purposes, which becomes computationally time exhaustive further affecting the learning procedure 61 . Therefore, the OSELM, which is an advanced form of ELM, operates in two learning stages utilizing the chunk-by-chunk approach i.e., initialization followed by the sequential learning stage. In the initialization phase, the H matrix is packed like ELM, for later usage. The randomized weights together with the biases are allocated to respective chunks of primary wheat yield (W) data to determine the output matrix of the www.nature.com/scientificreports/ OSELM hidden layers. Then the sequential learning stage is launched either in a one-by-one manner or a lumpby-lump fashion where the one-time data utilization is not permissible. More specific information on OSELM can be found in the following (e.g., [62][63][64]. For a given training set k−1 in the initialization phase: The term k−1 indicates the training dataset whereas x j is the input data and t j is the jth parameter. The first output weight is provided by the following equation: The term ρ k−1 is showing the initial output weight, t is denoting the hidden layer output matrix, and t k−1 = [t 1 , . . . , t k−1 ] t is the training data matrix. The biases and random weights are allocated in a small chunk in the initialization stage to calculate the hidden layer output matrix in the initial wheat yield (W) training data. The sequential learning phase is then commenced where RLS algorithm is utilized to modify the output weights in a recursive way 61 . The output weights in OSELM are recursively updated based on the intermediate results in the last iteration and the newly arrived data, which is deleted immediately once the features have been learnt, and therefore, the calculation overhead and the memory requirement of the algorithm are significantly decreased.

Random forest (RF).
The RF model is a regression tree-based learning approach whereby the bootstrapping and bagging are the underpinning modeling techniques on which the RF ensemble modeling approach is constructed upon 65,66 . Using a random bagging technique, the RF model develops ensembles in which each node is linked randomly by choosing the relevant inputs for better efficiency while avoiding overfitting 67 . The subsequent stages provide a brief account of the RF model designing: i. Perform bootstrapping on the input predictors to create n-bootstrapped based trees (i.e., n trees ). ii. Determine the largest no. of split variables (m try ) by means of random sampling with a non-prune regression tree. iii. Aggregate the simulated n trees to predict wheat yield (W).
The RF algorithm has been used in many applications including water quality 46 , soil moisture 68 , ecological 69 , hydrological 70 and solar radiation 71,72 forecasting applications.

Case study description and data
Study region and wheat yield data. Pakistan's Federal Bureau of Statistics and the Agriculture Marketing Information Services, Directorate of Agriculture provided the wheat yield data (Economics & Marketing) 73,74 . The study stations included are of high agricultural significance for the Punjab province, Pakistan. In this paper, each district is represented as a station. Agricultural sectors in Punjab province play an important part with economic contributions of 56.1-61.5% 75 . Further, a widespread irrigation network enables this province's rich agricultural zone. Considering Punjab as a key agronomic belt, the adoption of AI techniques to predict wheat yield is an interesting research endeavour. To establish the time series wheat yield dataset, the district-level wheat production was assimilated. To construct this dataset, the areal (district level) productions of wheat were acquired through the provincial Crop Reporting Services which had been compiled by the Economic Wing of the devolved Ministry of Food and Agriculture, and later by the Pakistan Federal Bureau of Statistics. The acquired data had some missing wheat yield values for 2009. To overcome this issue, the average of all other data for the period 1981-2013 is used to recover the missing data of the predictor and the corresponding target stations. Figure 1a,b show the provinces in Pakistan and the map of all districts in Punjab province (current study region). The figure illustrates the study stations (i.e., the major districts) of wheat farming. Figure 2a-f present the total of 6 maps which represents the testing stations (yellow colour), training stations (red colour), and the stations where wheat yield data is not available (green colour) and the stations that are not selected by ACO algorithm (blue colour). Figures 1 and 2 were generated using the GIS software 76 . A total of 27 stations were considered with data from 1981 to 2013. To obtain the wheat yield time series, out of 27 stations, 26 stations were used for the selection of the best stations for training to develop the model in relation to the remaining (1) testing station. Each time, 26 stations were used to select the best stations for training subsets against the 6 test stations. Table 1 presents basic statistics (i.e., latitude, longitude, and elevations), maximum, minimum, standard deviation, skewness, and kurtosis) of the present study stations.
Design of two-phase hybrid ACO-OSELM model. The two-phase hybrid ACO-OSELM model was developed on MATLAB R2016b platform, (The Math Works Inc. USA) with Pentium 4 dual-core Central Processing Unit (CPU). To develop the proposed two-phase universal ACO-OSELM model, historical wheat yield data series were used. In this paper, W represents the wheat yield, W obs denotes the rvobserved wheat yield while W pred represents the predicted wheat yield. The original wheat yield data that had statistically significant lagged values at (t − 1) were employed as the input predictors in the first phase of the developmental stages. The development of the two-phase hybrid ACO-OSELM model involved the following phases: www.nature.com/scientificreports/ Phase 1. In the foremost phase, the selected stations for model training were determined using the robust ACO feature selection strategy. Then, the user-defined parameters were defined with the ant numbers as 10 having 20 iterations and the initial pheromone factor was 1. For every station, the number of predictor stations (features) was defined prior to running the model. For Rahimyar Khan, the number of selected feature stations was 22, D. G. Khan (20), Kasur (19), Sialkot (17), Rawalpindi (12) and Jhang (14). The selected training stations with their correlation r against testing stations are described in Table 2. www.nature.com/scientificreports/ The proposed two-phase hybrid ACO-OSELM model was trained and tested on a longer time series as well as at a shorter time series to assess the accuracy and universal performance of the model. This is to ensure that the model could be applied to other locations in Pakistan. In addition, a different number of surrounding predictor stations (features) were defined for every other testing station. Essentially, the Rahimyar Khan testing station had the largest number of surrounding training stations (i.e., 22) having the longest time series with 726 data points in the time series. On the other hand, Rawalpindi testing station has 12 training stations being selected with 396 data points being the shortest time-series used in the study. Table 3 shows the training data lengths for respective stations. In addition, the pheromone exponential weights and heuristic exponential weights were kept as 1. Figure 3 plots the RMSE errors that occurred in optimizing the cost and objective function of the ACO algorithm for identifying the best feature stations.

Scientific Reports
Since the data consisted of annual values from 1981 to 2013, which resulted in a total of 33 data records. ML models sometimes perform poorly on shorter time series. To handle this issue, we adopted the approach of selection of stations by ACO algorithm for training purposes and test the proposed model on the whole time-series data for respective testing stations. This technique of utilizing yield data from surrounding study stations to predict the yield of test stations are practically useful since it can enable the modellers to extract similar features and patterns prevalent at a predictor station to be analysed to estimate the yield at a testing station. This modelling approach does not required the splitting of the data in the traditional manner.
After carefully selecting the training stations for respective testing stations using the ACO algorithm, their correlation r (of selected training stations) against testing stations were calculated to confirm the linear relationship among them. For the study station Rahimyar Khan, the training station Khanewal registered the highest value of r ≈ 0.855, followed by Bahawal Nagar (r ≈ 0.854). Similarly, Muzaffar Garh (r ≈ 0.881) and Rajanpur (r ≈ 0.861) have the largest values of correlation with station D. G. Khan. For the study station, Kasur, Gujranwala, and Shekhupura attained the highest values of (r ≈ 0.950, 0.947). Table 2 summarizes all the correlation coefficients for respective stations. On the other hand, Station Kasur has the smallest cost to objective function RMSE followed by Jhang station (Fig. 3). Table 3 presents the number of datum points for training and testing purposes in each station with a ratio of selected stations against testing stations, with training and testing data distribution parameters. To prevent the differences in skewness in training and testing affecting the outcomes, data normalization was carried out using the following equation: Table 1. Geographic properties and wheat yield statistics of the study stations for Punjab, Pakistan.

Geographic characteristics
Wheat yield statistics (kg ha −1 ) www.nature.com/scientificreports/  Table 3. Training data points (in terms of selected training stations) and testing data points for each testing station using the ACO algorithm. The ratio of selected stations against testing stations, skewness, and kurtosis of training and testing data. www.nature.com/scientificreports/ In Eq. (9) W indicates input/output of the wheat yield data, W min is the smallest value, W max is the largest value of wheat yield in the dataset, and W norm is the desired normalized value. The normalization process overcomes data fluctuations caused by inherent data features/patterns 77 . It essentially is invertible and in no way affects the results 77 . Figure 4 presents the time series of the tested study stations constructed from the selected features using the ACO algorithm.  www.nature.com/scientificreports/ Phase 2. The partial autocorrelation function (PACF) was employed to calculate and determine the statistically significant lags of historical wheat yield data series as in Fig. 5. Moreover, the cross-validation or any data randomized approach cannot be adopted as time-series data by definition occur in a temporal order/sequence and this order or sequence must be preserved in order to keep the structure of the series intact 78 . These significant lagged inputs at (t − 1) were incorporated as the input predictor in the OSELM model to forecast the yield W. Different activation functions (sigmoid, sine, hardlim, radial basis) were tested to determine the best activation function and the radial basis (rbf) and sigmoid (sig) functions were found to be the optimal ones. Consequently, different combinations of hidden neuron ranging from 7 to 35 were trialed with block size being fixed at 100. The second significant lag (t − 2) was also utilized in the proposed two-phase hybrid ACO-OSELM model to check whether model performance increases or not. But upon utilizing the lag (t − 2), the accuracy of the proposed two-phase hybrid ACO-OSELM model decreased, so it was not considered in this paper. Similarly, the benchmark models ELM and RF models were developed resulting in ACO-ELM and ACO-RF models respectively. Figure 6 displays the model schematics. Then model training performances of the proposed hybrid two-phase ACO-OSELM were assessed via correlation coefficient (r) and root-mean-squarederror (RMSE) as shown in Table 4.

Latitude (N) Longitude (E) Elevation (m) Mean
The r and RMSE values attained during the training period by the two-phase hybrid ACO-OSELM models for wheat yield prediction at Rahimyar Khan and D. G. Khan were seen to be: (r = 0.812, 0.790, RMSE = 374.82, 381.57 kg ha −1 ). Equivalent metrics for Kasur and Sialkot were found to be: (r = 0.804, 0.798, RMSE = 370.49, 386.18 kg ha −1 ) and finally for Rawalpindi and Jhang were: (r = 0.832, 0.799, RMSE = 356.80, 353.55 kg ha −1 ). In addition, the training performances of comparative ACO-ELM and ACO-RF models were also studied ( Table 4). The performance of the proposed two-phase hybrid ACO-OSELM model was relatively high during the training phase, and it is conjectured that the ACO-OSELM model accuracy in the testing phase for wheat yield prediction at these tested stations would be higher as well.
Setting and tuning parameter optimization. To attain optimum precision, one of the most crucial tasks in designing the prediction model is to adjust the tuning and pruning of parameters associated with the models. Various approaches are adopted to fine-tune the parameters to acquire the desired optimum model. The trial and error strategy was utilized to get the optimum parameters during the constructing phase of the ACO-OSELM, ACO-ELM, and ACO-RF model to predict the wheat yield 79 . Table 5   www.nature.com/scientificreports/ split predictors. The details on fine-tuning of these parameters are provided in Table 5 for optimally performing ACO-OSELM, ACO-ELM, and ACO-RF model for all selected study stations. iii. Nash-Sutcliffe coefficient (NS E ): iv. Root mean square error (RMSE, kg ha −1 ): Table 4. Training performance of two-phase hybrid ant colony optimization algorithm coupled online sequential extreme learning machines (ACO-OSELM) versus ACO-ELM and ACO-RF models with correlation coefficient (r) and root mean squared error (RMSE, kg ha −1 ).  Wpred,i represents respective observed and predicted averages of W while N is the number of data points in the testing phase. The value of correlation coefficient (r) can be in the range of − 1 and + 1, which demonstrates the associations in terms of the proportion of variance in between the observed and predicted W from the machine learning models 82 .

Testing stations Lags
A value of + 1 shows that the observed and forecasted values are highly correlated with the least variances. The second metric Willmott's Index (WI) ranges between 0 and 1. The WI overcomes the insensitivity issues as the differences between the observed and forecasted values are not squared and the ratio of the mean squared error in place of the differences are considered in computations 87,88 . The next metric, i.e., Nash-Sutcliffe Efficiency (NS E ) has the ideal value of 1 and spans till negative infinity that essentially compares the variance of predicted with the observed values 89 . In addition, the computation of error measures RMSE and MAE are based on the aggregation of residuals of observed and predicted W values 90 . The higher W values are largely captured by the RMSE whereas the MAE equally assesses all variations regardless of the sign, yet both range from 0 (ideal value) to positive infinity. On the other hand, the Legates-McCabe's (LM) index is a more robust tool developed to address the limitations of both the W and NS E and the value is bound between 0 and 1 (the ideal value) 91 .
However, these metrics should not essentially be used to compare model performance at geographically diverse stations 92 , as these metrics are in their absolute terms. As such the relative values of root mean square error (RRMSE) and mean absolute error (RMAE) were utilized for this purpose 93 . The relative values are in percentages and for a model to be rated as outstanding, the (RRMSE, RMAPE) < 10%. For models rated as good the range is 10% < (RRMSE, RMAE) < 20%, while the model is fair if 20% < (RRMSE, RMAE) < 30% and if the (RRMSE, RMAE) > 30% the model is considered to have poor prediction performance 94,95 .

Modelling results and analysis
The proposed two-phase hybrid ACO-OSELM is evaluated with ACO-ELM and ACO-RF models, based on evaluation metrics ("Setting and tuning parameter optimization" section), diagnostic plots together with error distributions. Figure 7 displays a scatterplot with the respective coefficients of determination (r 2 ) depicting the level of associations between the predicted and observed wheat yield (W) overlayed with the goodness-of-fit line and the linear equation. Essentially, the closer the linear equation is to the y = mx representation and the closer the r 2 is to unity, the better the model performance is. The proposed two-phase hybrid ACO-OSELM model has better predictive potential than ACO-ELM and ACO-RF models in terms of r 2 (ACO-OSELM ≈ 0.995, ACO-ELM ≈ 0.996, ACO-RF ≈ 0.862) for Kasur. Again, the proposed two-phase hybrid ACO-OSELM model is more accurate for Sialkot, r 2 (ACO-OSELM ≈ 0.974, ACO-ELM ≈ 0.936, ACO-RF ≈ 0.892), and Rawalpindi stations in terms of the achieved r 2 (ACO-OSELM ≈ 0.945, ACO-ELM ≈ 0.924, ACO-RF ≈ 0.814). The proposed twophase hybrid ACO-OSELM model for other stations Rahimyar Khan, D. G. Khan, and Jhang is reasonably good compared to ACO-ELM and ACO-RF models (Fig. 7). The better accuracy of the proposed two-phase hybrid ACO-OSELM model against the comparison models for all the study regions is confirmed by the linear regression equation and the goodness-of-fit in addition to attaining larger r 2 values.
Moreover, comparative boxplots of the proposed two-phase hybrid ACO-OSELM model with ACO-ELM and ACO-RF models for each station were established.  www.nature.com/scientificreports/ benchmark models. Similarly, the proposed two-phase hybrid ACO-OSELM model performed well for Sialkot and Kasur stations in predicting wheat yield followed by the ACO-ELM and ACO-RF models. The boxplot in Fig. 8 clearly shows that the proposed two-phase hybrid ACO-OSELM model at all six stations outperformed the comparative models. The preciseness of the proposed two-phase hybrid ACO-OSELM is further evaluated with comparative ACO-ELM and ACO-RF models based on r, RMSE, and MAE ( Table 6) can be considered a better data-intelligent tool for wheat yield prediction in contrast to the ACO-ELM and ACO-RF models. A vector field evaluation (VFE) diagram (Fig. 9) presents a more concise statistical summary of the associations of predicted and observed wheat yield matched based upon the respective WI values. A VFE diagram is a generalization of the Taylor diagram that provides an evaluation of model performances in terms of vector fields 96 . For Rahimyar Khan, the WI of the proposed two-phase hybrid ACO-OSELM model with observations was ~ 0.98, followed by ACO-ELM ≈ 0.97 and ACO-RF ≈ 0.87. The WI ~ 0.99 of the ACO-OSELM model was closest to the observed wheat yield as compared to ACO-ELM and ACO-RF for D. G. Khan stations. Similarly, the proposed two-phase hybrid ACO-OSELM models were found to be the best performing models for Kasur station (WI ≈ 0.98) that were within close proximity of the observed wheat yield followed by ACO-ELM (≈ 0.97) and ACO-RF (≈ 0.92) models. For other stations Sialkot and Jhang, the proposed two-phase hybrid ACO-OSELM model is closer to the observed W as compared to the ACO-ELM and ACO-RF models. For the Rawalpindi station, the ACO-RF model was marginally better than ACO-OSELM and ACO-ELM models. Overall, the WI of the proposed two-phase ACO-OSELM model was closely distributed to the observed baseline compared to the other models.  www.nature.com/scientificreports/  (Table 7) in comparison to the counterpart models revealing the better performances of the proposed two-phase hybrid ACO-OSELM models.
Furthermore, the empirical cumulative distribution function (ECDF, Fig. 10) at all stations depicts that the proposed two-phase hybrid ACO-OSELM method was reasonably better and superior to both the ACO-ELM and ACO-RF models. Based on the error (0 to ± 400 kg ha −1 ) for the Rahimyar Khan, D. G. Khan, and Kasur station, (0 to ± 600 kg ha −1 ) for Rawalpindi and Jhang station while (0 to ± 1000 kg ha −1 ) for Sialkot station, Fig. 10 clearly proves that the proposed two-phase hybrid ACO-OSELM method was the most precise model in predicting wheat yield.   Table 8. It demonstrated that D. G. Khan is the station where the ACO-OSELM wheat yield predicting model performed the best with RRMSE ≈ 3.00 and RMAE ≈ 2.25%. The ACO-OSELM model was found to generate the least relative percentage error values (i.e., RRMSE, RMAE) at all tested stations except for the Rawalpindi station. However, the predicted errors generated by the proposed two-phase hybrid ACO-OSELM model were low in terms of their relative error values, and more importantly, the error values were within the recommended range of 10% threshold for an excellent model classification, except for Rawalpindi station 97 . Figure 11 presents the absolute prediction error |PE| in each year from 1981-2013 of the proposed two-phase hybrid ACO-OSELM versus ACO-ELM and ACO-RF models at the testing stations in the form of polar plots. For all stations, the prediction errors generated by the proposed two-phase hybrid ACO-OSELM were very low compared to the ACO-ELM and ACO-RF models. This was justified by the minimum values of relative prediction errors. The |PE| errors were significantly smaller in each year for the proposed two-phase hybrid ACO-OSELM model as compared to ACO-ELM and ACO-RF models in Rahimyar Khan, D. G. Khan, Kasur, Sialkot, Rawalpindi and Jhang stations. Overall, the proposed two-phase hybrid ACO-OSELM model generated better significant accuracy with smaller error statistics and higher WI.

Discussions
Developing strategies for accurate crop yield prediction that can address food scarcity issues, decision-making on national imports and exports, and setting the prices in agriculture markets, can play an important role in policymaking, particularly in agriculture-based nations such as Pakistan. This study was aimed at designing a novel two-phase hybrid ACO-OSELM model using significant lag at (t − 1) to predict future wheat yield. This is a practically useful approach for crop yield management in terms of using the wheat yield data from several nearby stations in developing better agricultural practices and efficient precision agricultural technologies. For example, the methodology can be used in remote areas where meteorological data is not available due to limited resources. The research framework in this study can be applied to any other study station where yield data (whether it is wheat or any other crop) are available from surrounding stations to provide an accurate prediction.
The proposed two-phase hybrid ACO-OSELM model with its counterpart models (ACO-ELM and ACO-RF) was suitably evaluated, which revealed smaller relative percentage errors in terms of RRMSE and RMAE being generated. Respectively, reasonably large statistical correlation metric values of Legates-McCabe's between predicted and observed yielded for D. G. Khan and other tested stations (Tables 7, 8). The model performances were high since the percentage errors achieved were less than 10%. Thus, the proposed two-phase hybrid ACO-OSELM model can suitably be used to predict the wheat yield where the prediction of a crop commodity is likely to become more important due to ever-increasing demand.
The proposed two-phase hybrid OSELM model can assist the government's national policymakers and agricultural engineers in minimizing uncertainties in crop estimates reducing price hikes as well as unwarranted wastages 98 . Since the proposed two-phase hybrid ACO-OSELM model offers better forecasting potential together with being fast and robust, it can possibly be explored in predicting other crop yields including Rice, Maize, Cotton, Sugarcane, and Oilseeds to generate similar optimal predictions in follow-up studies.
The utilization of historical wheat yield data as inputs to predict the future yield carries some limitations. Certainly, weather conditions are a big driver for any agricultural yield. Hence, to enhance the scope of future studies, predictor inputs consisting of meteorological data (i.e., precipitation, air temperatures, soil moisture, wind speed, solar radiation, etc.) need to be used to predict future crop yield as crop production amounts are largely contingent upon these parameters. Satellite-based remotely sensed data and/or data from atmospheric simulation models (e.g., 31,99,100 ) as predictor variables are also likely to add great value to the crop yield modeling in remote agricultural areas with limited datasets. Incorporation of remotely sensed photosynthetically active radiation (PAR) that reveals crop health could also be focussed on, in an independent study. In addition, fertilizer/manure usage data accompanied by relevant soil characteristics (e.g., texture, pedality, hydraulics, porosity, bulk density, thickness, and soil organic carbon content) could also be explored in the proposed twophase hybrid ACO-OSELM model. Fields that use irrigation for production need to utilize irrigation statistics to improve crop yield predictions. Further, as process-based modeling is resource-demanding which emerging Table 8. Geographic comparison of the accuracy of the ACO-OSELM versus ACO-ELM and ACO-RF models in terms of relative root mean squared error (RRMSE, %) and the relative mean absolute error (RMAE, %) computed within the test stations. Note that the best model is boldfaced (underline). www.nature.com/scientificreports/ www.nature.com/scientificreports/ agricultural nations like Pakistan unable to afford, the proposed two-phase hybrid ACO-OSELM model could be used as a feasible option. An ensemble modeling approach could further improve the two-phase hybrid ACO-OSELM modeling with the possibility of achieving better results. Ensemble modeling would provide better confidence in predictions making the model more reliable in strategic decision-making, as uncertainties between several forecasted data would be captured and displayed in the outputs. Quantum-Behaved PSO and the Firefly Algorithms could also be used to identify training stations that have been tested to hybridize with the OSELM (e.g., 101,102 ). Future works could apply, empirical wavelet transform-EWT 103 , empirical mode decomposition-EMD 104 , and singular value decomposition-SVD 105 , as data pre-processing tools in modeling and predicting crop yields.

Conclusions
The current study was adopted to develop a robust two-phase hybrid ACO-OSELM model to predict wheat yield. Lagged wheat yields from several neighbouring stations were used for training purposes as model predictors for the candidate station. Wheat yield data for the period of 1981-2013 from 26 stations were pooled and the best training stations were selected by the ACO algorithm on the basis of feature selection corresponding to the 27th test station. The selected feature stations were used to construct a time series where the significant lags at (t − 1) were used to develop the proposed two-phase hybrid ACO-OSELM model to achieve better accuracy. Several evaluation criteria including diagnostic plots were adopted to judge the accuracy of the proposed two-phase hybrid ACO-OSELM model. The proposed hybrid ACO-OSELM outperformed the counterpart models for wheat yield prediction. The prediction errors metrics for the best station D. G. Khan registered by ACO-OSELM model were RMSE ≈ 67.12 kg ha −1 and MAE ≈ 42.37 kg/ha. The normalized performance metrics for the D. G. Khan station (r ≈ 0.997, WI ≈ 0.989 and NS E ≈ 0.978. The performance assessment of the ACO-OSELM model in relation to ACO-ELM and the ACO-RF models via Legates-McCabe's indices were in agreement. The LM values between the predicted and observed wheat yield for the D. G. Khan study station were LM ≈ 0.884 (ACO-OSELM), 0.879 (ACO-ELM) and 0.612 (ACO-RF), respectively and the relative errors, RRMSE and RMAE were very small: 3.00%, 2.25% (ACO-OSELM) compared with 3.09%, 2.39% (ACO-ELM) and 8.48%, 7.33% (ACO-RF). Since the relative percentage errors, RRMSE and RMAE showed that at D. G. Khan station the proposed two-phase hybrid ACO-OSELM model performed the best as compared to other stations, evidently geographic variability does influence the outcomes to a certain degree. This essentially is a baseline study whereby wheat yield data from several stations are being utilized to predict wheat yield more accurately that can potentially be extended to forecasting using other climatological parameters in future studies. Similarly, other agricultural crop yield predictions could be explored with the proposed two-phase hybrid ACO-OSELM model that will assist policymakers and decision-makers in the better management of crop yield and price predictions. More importantly, accurate wheat and other crop yield predictions can be used to alert impacted stakeholders and the government to avert food security issues.